Part 4

RNA velocity

library(velocyto.R)
library(SeuratWrappers)
library(Seurat)
source("tianfengRwrappers.R")
preprocess <- function(loom_path, n){
  velodata <- ReadVelocity(file = loom_path)
  
  #匹配两次的barcode
  tfunc <- function(s)
  {
    s <- strsplit(s,".*:",fixed = F)[[1]][2]
    s <- strsplit(s,"x",fixed = T)[[1]]
    s <- paste(s,as.character(n),sep = "-")
    return(s)
  }
  
  velodata[["spliced"]]@Dimnames[[2]] = as.character(lapply(velodata[["spliced"]]@Dimnames[[2]],tfunc))
  velodata[["unspliced"]]@Dimnames[[2]] = as.character(lapply(velodata[["unspliced"]]@Dimnames[[2]],tfunc))
  velodata[["ambiguous"]]@Dimnames[[2]] = as.character(lapply(velodata[["ambiguous"]]@Dimnames[[2]],tfunc))
  return(velodata)
}
ctrl <- readRDS("ctrl.rds")
ctrl_velo <- as.Seurat(x = preprocess("./HSD0524.loom",1))
ctrl_velo <- subset(ctrl_velo, cells = WhichCells(ctrl))

ctrl_velo <- ctrl_velo %>% 
    PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F,assay = "spliced") %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)


# for(i in levels(Idents(ctrl))){
#   (DimPlot(ctrl_velo, reduction = "umap",cells.highlight = WhichCells(ctrl, expression = `celltype` == i), pt.size = 0.5, label = T)+ggtitle(i)) %>% print()
# }
ctrl_velo@reductions[["umap"]] <- ctrl@reductions[["umap"]]
ctrl_velo@reductions[["umap"]]@cell.embeddings <-
  ctrl@reductions[["umap"]]@cell.embeddings[colnames(ctrl_velo),]

DimPlot(ctrl_velo, reduction = "umap",group.by = 'seurat_clusters', pt.size = 0.5, label = T)

ctrl_velo$Classification1 <- ctrl$celltype
Idents(ctrl_velo) <- ctrl_velo$Classification1
ctrl_velo@assays[["RNA"]] <- ctrl@assays[["SCT"]]

# saveRDS(ctrl_velo,"ctrl_velo.rds") 
# ctrl_velo$Classification1 <- Idents(ctrl)
ident.colors <- colors_list[1:length(x = levels(x = ctrl_velo))]
names(x = ident.colors) <- levels(x = ctrl_velo)
cell.colors <- ident.colors[Idents(object = ctrl_velo)]

names(x = cell.colors) <- colnames(x = ctrl_velo)

png("velocity.png")
show.velocity.on.embedding.cor(emb = Embeddings(object = ctrl_velo, reduction = "umap"), vel = Tool(object = ctrl_velo, slot = "RunVelocity"), n = 200, scale = "sqrt", cell.colors = ac(x = cell.colors, alpha = 0.5), 
    cex = 0.8, arrow.scale = 3, show.grid.flow = TRUE, min.grid.cell.mass = 0.5, grid.n = 40, arrow.lwd = 1, 
    do.par = FALSE, cell.border.alpha = 0.1)
dev.off()
h5ad_output <- "hsd_velo.h5ad"  

library(reticulate)
library(Matrix)
writeMM(t(ctrl_velo@assays$RNA@counts), file='combined.mtx')
writeMM(t(ctrl_velo@assays$spliced@counts), file='spliced.mtx')
writeMM(t(ctrl_velo@assays$unspliced@counts), file='unspliced.mtx')
write.csv(rownames(ctrl_velo@assays$spliced@counts), file = "genes.csv", row.names = FALSE)
write.csv(ctrl_velo@reductions$umap@cell.embeddings, file = "umap.csv", row.names = FALSE)
write.csv(ctrl_velo@reductions$pca@cell.embeddings, file = "pca.csv", row.names = FALSE)

write.csv(colnames(ctrl_velo@assays$spliced@counts), file = "cells.csv", row.names = FALSE)
write.csv(ctrl_velo@meta.data, file = "meta.csv", row.names = FALSE)
source_python('~/scRNAseq/vascular-analysis/build.py')
build(h5ad_output, pca = TRUE, umap = TRUE)
filename = "hsd_velo.h5ad"  
import anndata
import pandas as pd
from itertools import chain
#import scvelo as scv

ad = anndata.read_mtx('spliced.mtx')
ad.layers['spliced'] = anndata.read_mtx('spliced.mtx').X
ad.layers['unspliced'] = anndata.read_mtx('unspliced.mtx').X
ad.obs = pd.read_csv('meta.csv',header=0)
ad.obsm['X_umap'] = pd.read_csv('umap.csv',header=0).values
ad.obsm['X_pca'] = pd.read_csv('pca.csv',header=0).values


ad.var_names = pd.Index(list(chain.from_iterable(pd.read_csv('genes.csv', header=0,dtype='string').values.tolist()))).astype(str)
ad.obs_names = list(chain.from_iterable(pd.read_csv('cells.csv', header=0).values.tolist()))

ad.write(filename)
file.remove('combined.mtx')
[1] TRUE
file.remove('spliced.mtx')
[1] TRUE
file.remove('unspliced.mtx')
[1] TRUE
file.remove('genes.csv')
[1] TRUE
file.remove('cells.csv')
[1] TRUE
file.remove('umap.csv')
[1] TRUE
file.remove('pca.csv')
[1] TRUE
file.remove('meta.csv')
[1] TRUE

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

LS0tCnRpdGxlOiAiUk5BX3ZlbG9jaXR5IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIFBhcnQgNAojIyBSTkEgdmVsb2NpdHkKCmBgYHtyfQpsaWJyYXJ5KHZlbG9jeXRvLlIpCmxpYnJhcnkoU2V1cmF0V3JhcHBlcnMpCmxpYnJhcnkoU2V1cmF0KQpzb3VyY2UoInRpYW5mZW5nUndyYXBwZXJzLlIiKQpgYGAKYGBge3J9CnByZXByb2Nlc3MgPC0gZnVuY3Rpb24obG9vbV9wYXRoLCBuKXsKICB2ZWxvZGF0YSA8LSBSZWFkVmVsb2NpdHkoZmlsZSA9IGxvb21fcGF0aCkKICAKICAj5Yy56YWN5Lik5qyh55qEYmFyY29kZQogIHRmdW5jIDwtIGZ1bmN0aW9uKHMpCiAgewogICAgcyA8LSBzdHJzcGxpdChzLCIuKjoiLGZpeGVkID0gRilbWzFdXVsyXQogICAgcyA8LSBzdHJzcGxpdChzLCJ4IixmaXhlZCA9IFQpW1sxXV0KICAgIHMgPC0gcGFzdGUocyxhcy5jaGFyYWN0ZXIobiksc2VwID0gIi0iKQogICAgcmV0dXJuKHMpCiAgfQogIAogIHZlbG9kYXRhW1sic3BsaWNlZCJdXUBEaW1uYW1lc1tbMl1dID0gYXMuY2hhcmFjdGVyKGxhcHBseSh2ZWxvZGF0YVtbInNwbGljZWQiXV1ARGltbmFtZXNbWzJdXSx0ZnVuYykpCiAgdmVsb2RhdGFbWyJ1bnNwbGljZWQiXV1ARGltbmFtZXNbWzJdXSA9IGFzLmNoYXJhY3RlcihsYXBwbHkodmVsb2RhdGFbWyJ1bnNwbGljZWQiXV1ARGltbmFtZXNbWzJdXSx0ZnVuYykpCiAgdmVsb2RhdGFbWyJhbWJpZ3VvdXMiXV1ARGltbmFtZXNbWzJdXSA9IGFzLmNoYXJhY3RlcihsYXBwbHkodmVsb2RhdGFbWyJhbWJpZ3VvdXMiXV1ARGltbmFtZXNbWzJdXSx0ZnVuYykpCiAgcmV0dXJuKHZlbG9kYXRhKQp9CmBgYAoKYGBge3J9CmN0cmwgPC0gcmVhZFJEUygiY3RybC5yZHMiKQpjdHJsX3ZlbG8gPC0gYXMuU2V1cmF0KHggPSBwcmVwcm9jZXNzKCIuL0hTRDA1MjQubG9vbSIsMSkpCmN0cmxfdmVsbyA8LSBzdWJzZXQoY3RybF92ZWxvLCBjZWxscyA9IFdoaWNoQ2VsbHMoY3RybCkpCgpjdHJsX3ZlbG8gPC0gY3RybF92ZWxvICU+JSAKICAgIFBlcmNlbnRhZ2VGZWF0dXJlU2V0KHBhdHRlcm4gPSAiXk1ULSIsIGNvbC5uYW1lID0gInBlcmNlbnQubXQiKSAlPiUKICAgIFNDVHJhbnNmb3JtKHZhcnMudG8ucmVncmVzcyA9ICJwZXJjZW50Lm10IiwgdmVyYm9zZSA9IEYsYXNzYXkgPSAic3BsaWNlZCIpICU+JSAKICAgIFJ1blBDQSgpICU+JSBGaW5kTmVpZ2hib3JzKGRpbXMgPSAxOjIwKSAlPiUgCiAgICBSdW5VTUFQKGRpbXMgPSAxOjIwKSAlPiUgCiAgICBGaW5kQ2x1c3RlcnMocmVzb2x1dGlvbiA9IDAuMSkKCgojIGZvcihpIGluIGxldmVscyhJZGVudHMoY3RybCkpKXsKIyAgIChEaW1QbG90KGN0cmxfdmVsbywgcmVkdWN0aW9uID0gInVtYXAiLGNlbGxzLmhpZ2hsaWdodCA9IFdoaWNoQ2VsbHMoY3RybCwgZXhwcmVzc2lvbiA9IGBjZWxsdHlwZWAgPT0gaSksIHB0LnNpemUgPSAwLjUsIGxhYmVsID0gVCkrZ2d0aXRsZShpKSkgJT4lIHByaW50KCkKIyB9CmBgYAoKCmBgYHtyfQpjdHJsX3ZlbG9AcmVkdWN0aW9uc1tbInVtYXAiXV0gPC0gY3RybEByZWR1Y3Rpb25zW1sidW1hcCJdXQpjdHJsX3ZlbG9AcmVkdWN0aW9uc1tbInVtYXAiXV1AY2VsbC5lbWJlZGRpbmdzIDwtCiAgY3RybEByZWR1Y3Rpb25zW1sidW1hcCJdXUBjZWxsLmVtYmVkZGluZ3NbY29sbmFtZXMoY3RybF92ZWxvKSxdCgpEaW1QbG90KGN0cmxfdmVsbywgcmVkdWN0aW9uID0gInVtYXAiLGdyb3VwLmJ5ID0gJ3NldXJhdF9jbHVzdGVycycsIHB0LnNpemUgPSAwLjUsIGxhYmVsID0gVCkKY3RybF92ZWxvJENsYXNzaWZpY2F0aW9uMSA8LSBjdHJsJGNlbGx0eXBlCklkZW50cyhjdHJsX3ZlbG8pIDwtIGN0cmxfdmVsbyRDbGFzc2lmaWNhdGlvbjEKY3RybF92ZWxvQGFzc2F5c1tbIlJOQSJdXSA8LSBjdHJsQGFzc2F5c1tbIlNDVCJdXQoKIyBzYXZlUkRTKGN0cmxfdmVsbywiY3RybF92ZWxvLnJkcyIpIAojIGN0cmxfdmVsbyRDbGFzc2lmaWNhdGlvbjEgPC0gSWRlbnRzKGN0cmwpCmBgYAoKCmBgYHtyfQppZGVudC5jb2xvcnMgPC0gY29sb3JzX2xpc3RbMTpsZW5ndGgoeCA9IGxldmVscyh4ID0gY3RybF92ZWxvKSldCm5hbWVzKHggPSBpZGVudC5jb2xvcnMpIDwtIGxldmVscyh4ID0gY3RybF92ZWxvKQpjZWxsLmNvbG9ycyA8LSBpZGVudC5jb2xvcnNbSWRlbnRzKG9iamVjdCA9IGN0cmxfdmVsbyldCgpuYW1lcyh4ID0gY2VsbC5jb2xvcnMpIDwtIGNvbG5hbWVzKHggPSBjdHJsX3ZlbG8pCgpwbmcoInZlbG9jaXR5LnBuZyIpCnNob3cudmVsb2NpdHkub24uZW1iZWRkaW5nLmNvcihlbWIgPSBFbWJlZGRpbmdzKG9iamVjdCA9IGN0cmxfdmVsbywgcmVkdWN0aW9uID0gInVtYXAiKSwgdmVsID0gVG9vbChvYmplY3QgPSBjdHJsX3ZlbG8sIHNsb3QgPSAiUnVuVmVsb2NpdHkiKSwgbiA9IDIwMCwgc2NhbGUgPSAic3FydCIsIGNlbGwuY29sb3JzID0gYWMoeCA9IGNlbGwuY29sb3JzLCBhbHBoYSA9IDAuNSksIAogICAgY2V4ID0gMC44LCBhcnJvdy5zY2FsZSA9IDMsIHNob3cuZ3JpZC5mbG93ID0gVFJVRSwgbWluLmdyaWQuY2VsbC5tYXNzID0gMC41LCBncmlkLm4gPSA0MCwgYXJyb3cubHdkID0gMSwgCiAgICBkby5wYXIgPSBGQUxTRSwgY2VsbC5ib3JkZXIuYWxwaGEgPSAwLjEpCmRldi5vZmYoKQoKYGBgCgpgYGB7cn0KaDVhZF9vdXRwdXQgPC0gImhzZF92ZWxvLmg1YWQiICAKCmxpYnJhcnkocmV0aWN1bGF0ZSkKbGlicmFyeShNYXRyaXgpCndyaXRlTU0odChjdHJsX3ZlbG9AYXNzYXlzJFJOQUBjb3VudHMpLCBmaWxlPSdjb21iaW5lZC5tdHgnKQp3cml0ZU1NKHQoY3RybF92ZWxvQGFzc2F5cyRzcGxpY2VkQGNvdW50cyksIGZpbGU9J3NwbGljZWQubXR4JykKd3JpdGVNTSh0KGN0cmxfdmVsb0Bhc3NheXMkdW5zcGxpY2VkQGNvdW50cyksIGZpbGU9J3Vuc3BsaWNlZC5tdHgnKQp3cml0ZS5jc3Yocm93bmFtZXMoY3RybF92ZWxvQGFzc2F5cyRzcGxpY2VkQGNvdW50cyksIGZpbGUgPSAiZ2VuZXMuY3N2Iiwgcm93Lm5hbWVzID0gRkFMU0UpCndyaXRlLmNzdihjdHJsX3ZlbG9AcmVkdWN0aW9ucyR1bWFwQGNlbGwuZW1iZWRkaW5ncywgZmlsZSA9ICJ1bWFwLmNzdiIsIHJvdy5uYW1lcyA9IEZBTFNFKQp3cml0ZS5jc3YoY3RybF92ZWxvQHJlZHVjdGlvbnMkcGNhQGNlbGwuZW1iZWRkaW5ncywgZmlsZSA9ICJwY2EuY3N2Iiwgcm93Lm5hbWVzID0gRkFMU0UpCgp3cml0ZS5jc3YoY29sbmFtZXMoY3RybF92ZWxvQGFzc2F5cyRzcGxpY2VkQGNvdW50cyksIGZpbGUgPSAiY2VsbHMuY3N2Iiwgcm93Lm5hbWVzID0gRkFMU0UpCndyaXRlLmNzdihjdHJsX3ZlbG9AbWV0YS5kYXRhLCBmaWxlID0gIm1ldGEuY3N2Iiwgcm93Lm5hbWVzID0gRkFMU0UpCnNvdXJjZV9weXRob24oJ34vc2NSTkFzZXEvdmFzY3VsYXItYW5hbHlzaXMvYnVpbGQucHknKQpidWlsZChoNWFkX291dHB1dCwgcGNhID0gVFJVRSwgdW1hcCA9IFRSVUUpCgpgYGAKCmBgYHtweXRob259CmZpbGVuYW1lID0gImhzZF92ZWxvLmg1YWQiICAKaW1wb3J0IGFubmRhdGEKaW1wb3J0IHBhbmRhcyBhcyBwZApmcm9tIGl0ZXJ0b29scyBpbXBvcnQgY2hhaW4KI2ltcG9ydCBzY3ZlbG8gYXMgc2N2CgphZCA9IGFubmRhdGEucmVhZF9tdHgoJ3NwbGljZWQubXR4JykKYWQubGF5ZXJzWydzcGxpY2VkJ10gPSBhbm5kYXRhLnJlYWRfbXR4KCdzcGxpY2VkLm10eCcpLlgKYWQubGF5ZXJzWyd1bnNwbGljZWQnXSA9IGFubmRhdGEucmVhZF9tdHgoJ3Vuc3BsaWNlZC5tdHgnKS5YCmFkLm9icyA9IHBkLnJlYWRfY3N2KCdtZXRhLmNzdicsaGVhZGVyPTApCmFkLm9ic21bJ1hfdW1hcCddID0gcGQucmVhZF9jc3YoJ3VtYXAuY3N2JyxoZWFkZXI9MCkudmFsdWVzCmFkLm9ic21bJ1hfcGNhJ10gPSBwZC5yZWFkX2NzdigncGNhLmNzdicsaGVhZGVyPTApLnZhbHVlcwoKCmFkLnZhcl9uYW1lcyA9IHBkLkluZGV4KGxpc3QoY2hhaW4uZnJvbV9pdGVyYWJsZShwZC5yZWFkX2NzdignZ2VuZXMuY3N2JywgaGVhZGVyPTAsZHR5cGU9J3N0cmluZycpLnZhbHVlcy50b2xpc3QoKSkpKS5hc3R5cGUoc3RyKQphZC5vYnNfbmFtZXMgPSBsaXN0KGNoYWluLmZyb21faXRlcmFibGUocGQucmVhZF9jc3YoJ2NlbGxzLmNzdicsIGhlYWRlcj0wKS52YWx1ZXMudG9saXN0KCkpKQoKYWQud3JpdGUoZmlsZW5hbWUpCmBgYApgYGB7cn0KZmlsZS5yZW1vdmUoJ2NvbWJpbmVkLm10eCcpCmZpbGUucmVtb3ZlKCdzcGxpY2VkLm10eCcpCmZpbGUucmVtb3ZlKCd1bnNwbGljZWQubXR4JykKZmlsZS5yZW1vdmUoJ2dlbmVzLmNzdicpCmZpbGUucmVtb3ZlKCdjZWxscy5jc3YnKQpmaWxlLnJlbW92ZSgndW1hcC5jc3YnKQpmaWxlLnJlbW92ZSgncGNhLmNzdicpCmZpbGUucmVtb3ZlKCdtZXRhLmNzdicpCmBgYAoKQWRkIGEgbmV3IGNodW5rIGJ5IGNsaWNraW5nIHRoZSAqSW5zZXJ0IENodW5rKiBidXR0b24gb24gdGhlIHRvb2xiYXIgb3IgYnkgcHJlc3NpbmcgKkN0cmwrQWx0K0kqLgoKV2hlbiB5b3Ugc2F2ZSB0aGUgbm90ZWJvb2ssIGFuIEhUTUwgZmlsZSBjb250YWluaW5nIHRoZSBjb2RlIGFuZCBvdXRwdXQgd2lsbCBiZSBzYXZlZCBhbG9uZ3NpZGUgaXQgKGNsaWNrIHRoZSAqUHJldmlldyogYnV0dG9uIG9yIHByZXNzICpDdHJsK1NoaWZ0K0sqIHRvIHByZXZpZXcgdGhlIEhUTUwgZmlsZSkuCgpUaGUgcHJldmlldyBzaG93cyB5b3UgYSByZW5kZXJlZCBIVE1MIGNvcHkgb2YgdGhlIGNvbnRlbnRzIG9mIHRoZSBlZGl0b3IuIENvbnNlcXVlbnRseSwgdW5saWtlICpLbml0KiwgKlByZXZpZXcqIGRvZXMgbm90IHJ1biBhbnkgUiBjb2RlIGNodW5rcy4gSW5zdGVhZCwgdGhlIG91dHB1dCBvZiB0aGUgY2h1bmsgd2hlbiBpdCB3YXMgbGFzdCBydW4gaW4gdGhlIGVkaXRvciBpcyBkaXNwbGF5ZWQuCg==